/*
 * Common routines used in the implementation
 * Salem Machaka
 * 04/03/98
 */

#include"mle.h"

/*
 * Function that returns Minimum
 */
double Min(double x1, double x2) {

    if (x1 < x2) return x1;
    else return x2;
}

/*
 * Indicator used in probelm formulation
 */
int Indicator(double x1, double M) {
    if (x1 <= (double) M) return (1);
    else return (0);
}

/*
 * Allocates memory for an integer vector from [beg...end]
 */
int *Int_Vector(int beg, int end) {
    int *vec;


    vec = (int *) malloc((end - beg + 1) * sizeof (int));
    if (vec == NULL) Abort("Error allocating memory in Int_vector, Aborting ...\n");
    return (vec - beg);

}

/*
 * Allocates memory for an double vector from [beg...end]
 */
double *Dec_Vector(int beg, int end) {
    double *vec;

    vec = (double *) malloc((end - beg + 1) * sizeof (double));
    if (vec == NULL) Abort("Error allocating memory in Dec_vector, Aborting ...\n");
    return (vec - beg);

}

/*
 * Free an int vector 
 */
void Free_Int(int *vec, int beg) {
    free(vec + beg);
}

/*
 * Free a dec vector 
 */
void Free_Dec(double *vec, int beg) {
    free(vec + beg);
}

/*
 * Allocate memory for an integer matrix [begrow...endrow][begcol...endcol]
 */
int **Int_Matrix(int begrow, int endrow, int begcol, int endcol) {

    int **mat;
    int i;


    mat = (int **) malloc((endrow - begrow + 1) * sizeof (int*));
    if (mat == NULL) Abort("Error allocating memory in Int_Matrix, Aborting ...\n");
    mat -= begrow;
    for (i = begrow; i <= endrow; i++) {
        mat[i] = (int *) malloc((endcol - begcol + 1) * sizeof (int));
        if (mat[i] == NULL) Abort("Error allocating memory in Int_Matrix, Aborting ...\n");
        mat[i] -= begcol;
    }
    return (mat);

}

/*
 * Allocate memory for a double matrix [begrow...endrow][begcol...endcol]
 */
double **Dec_Matrix(int begrow, int endrow, int begcol, int endcol) {
    double **mat;
    int i;


    mat = (double **) malloc((endrow - begrow + 1) * sizeof (double*));
    if (mat == NULL) Abort("Error allocating memory in Dec_Matrix, Aborting ...\n");
    mat -= begrow;
    for (i = begrow; i <= endrow; i++) {
        mat[i] = (double *) malloc((endcol - begcol + 1) * sizeof (double));
        if (mat[i] == NULL) Abort("Error allocating memory in Dec_Matrix, Aborting ...\n");
        mat[i] -= begcol;
    }
    return (mat);

}

/*
 * Free an int matrix
 */
void Free_Intmatrix(int **mat, int begrow, int endrow, int begcol) {
    int i;


    for (i = begrow; i <= endrow; i++)
        free(mat[i] + begcol);
    free(mat + begrow);

}




/*
 * Free a Dec matrix
 
void Free_Decmatrix(double **mat, int begrow, int endrow, int begcol, int endcol)
{
  int 		i;


  for(i=begrow; i<=endrow; i++)
        free(mat[i]+begcol);
  free(mat+begrow);

}
 */

/*
 * Abort in case of error with reporting the problem
 */
void Abort(char *message) {

    fprintf(stderr, "%s", message);
    exit(1);

}

/*
 * read in input files into Data Structures
 */
void Read_Files(int *nprobe, int *nclone, int *M, int *N,
        int ***Hyb, int **piorder, double *eata,
        double *rho, int *P, int argc, char **argv) {

    FILE *inp_param, *inp_conf, *inp_order, *inp_matrix;
    char mess[500];
    char filenames[5][100];
    int i, j;


    system("clear");
    /*
     * Checking input paramaters and existence of input files
     */
    if (argc != 2)
        Abort("Usage:\n mle <parameter_file> \n");




    inp_param = fopen(argv[1], "r");
    if (inp_param == NULL) {
        sprintf(mess, "Error opening %s. Aborting ...\n", argv[1]);
        Abort(mess);
    }

    /*
     * Reading parameters from parameter file
     */
    fscanf(inp_param, "%lf", &INIT_TEMP);
    fscanf(inp_param, "%d", &ITER);
    fscanf(inp_param, "%d", &SUCC);
    fscanf(inp_param, "%lf", &FACTOR);
    fscanf(inp_param, "%lf", &CRITERIA_CONJ);
    fscanf(inp_param, "%lf", &CRITERIA_BISECTION);
    for (i = 0; i < 5; i++)
        fscanf(inp_param, "%s", filenames[i]);
    fscanf(inp_param, "%d", &NUMB_ITER_LOG);
    fscanf(inp_param, "%d", &SEED);


    inp_conf = fopen(filenames[0], "r");
    if (inp_conf == NULL) {
        sprintf(mess, "Error opening %s. Aborting ...\n", filenames[0]);
        Abort(mess);
    }
    inp_matrix = fopen(filenames[1], "r");
    if (inp_matrix == NULL) {
        sprintf(mess, "Error opening %s. Aborting ...\n", filenames[1]);
        Abort(mess);
    }
    inp_order = fopen(filenames[2], "r");
    if (inp_order == NULL) {
        sprintf(mess, "Error opening %s. Aborting ...\n", filenames[2]);
        Abort(mess);
    }
    /*
     * Opening output & log files and Checking for error
     */
    OUT = fopen(filenames[3], "w");
    LOG = fopen(filenames[4], "w");




    /*
     * Reading variables from configuration file
     */
    fscanf(inp_conf, "%d", nprobe);
    fscanf(inp_conf, "%d", nclone);
    fscanf(inp_conf, "%d", M);
    fscanf(inp_conf, "%d", N);
    fscanf(inp_conf, "%lf", eata);
    fscanf(inp_conf, "%lf", rho);

    /*
     * Allocating memory and reading initial ordering (random) of probes
     */
    *piorder = Int_Vector(0, (*nprobe) + 1);
    for (j = 1; j <= *nprobe; j++)
        fscanf(inp_order, "%d", &(*piorder)[j]);
    (*piorder)[0] = (*piorder)[(*nprobe) + 1] = 0;


    /*
     * Allocating memory and reading hybridization data
     */
    *P = 0;
    *Hyb = Int_Matrix(1, *nclone, 1, *nprobe);
    for (i = 1; i <= *nclone; i++)
        for (j = 1; j <= *nprobe; j++) {
            fscanf(inp_matrix, "%1d", &(*Hyb)[i][j]);
            if ((*Hyb)[i][j] == 1) (*P)++;
        }

    /*
     * Closing the input files
     */
    fclose(inp_param);
    fclose(inp_conf);
    fclose(inp_order);
    fclose(inp_matrix);

}

void Initial_Param(int ***joint, double ***aa, int *gap, double *const1,
        double **y, int nprobe, int nclone, int N, int M,
        int **Hyb, double eata, double rho, int P) {
    int i, j, l;


    *y = Dec_Vector(1, nprobe + 1);
    /*
     * Initializing joint: joint[i][j] = # of clones that
     * hybridize to both i and j, joint[i][i] = 0, for all i
     */
    *joint = Int_Matrix(1, nprobe, 1, nprobe);
    for (i = 1; i <= nprobe; i++)
        for (j = 1; j <= nprobe; j++)
            (*joint)[i][j] = 0;
    for (i = 1; i <= nprobe - 1; i++)
        for (j = i + 1; j <= nprobe; j++) {
            for (l = 1; l <= nclone; l++)
                if ((Hyb[l][i] == 1) && (Hyb[l][j] == 1)) (*joint)[i][j]++;
            (*joint)[j][i] = (*joint)[i][j];
        }

    /*
     * Initializing aa: used in calculation
     */
    *aa = Dec_Matrix(1, nclone, 0, nprobe);
    for (i = 1; i <= nclone; i++)
        for (j = 0; j <= nprobe; j++)
            if (j == 0) (*aa)[i][j] = 0.0;
            else if (Hyb[i][j] == 0) (*aa)[i][j] = eata / (1 - rho);
            else (*aa)[i][j] = (1 - eata) / rho;

    *const1 = nclone * log(N - M) - P * log(rho / (1 - rho)) - nprobe * nclone * log(1 - rho);
    *gap = N - (nprobe * M);

}

/*
 * Function to evaluate the function value at point y
 */
double Func(double const1, double *R, double **aa, double *y, int M, int nprobe,
        int nclone, int *order) {

    /*
      double 	ret, sum;
      int 		i,j;
      char MYOUT[FILENAME_MAX] = "C:\\ONE.TXT";
      FILE *myout;

      myout = fopen(MYOUT,"w");
      ret = 0;
      for(i=1;i<=nclone;i++)
       {
         sum = 0.0;
         for(j=1;j<=nprobe+1;j++){
                       sum -= (1-aa[i][order[j]])*(1-aa[i][order[j-1]])*Min(y[j],(double)M);
                     }
         sum += R[i];
         ret -= log(sum);
         fprintf(myout,"R[i] %lf sum %lf ret %lf\n",R[i],sum,ret);
     
       }
       fclose(myout);
      ret += const1;



      return ret;
     */
    double ret, *sum;
    int i, j;
    //char MYOUT[FILENAME_MAX] = "C:\\ONE.TXT";
    //FILE *myout;

    //myout = fopen(MYOUT,"w");

    sum = Dec_Vector(1, nclone);

    for (i = 1; i <= nclone; i++) {
        sum[i] = 0.0;
        for (j = 1; j <= nprobe; j++) {
            sum[i] -= (1 - aa[i][order[j]])*(1 - aa[i][order[j - 1]]) * Min(y[j], (double) M);

        }
    }
    ret = 0;
    for (i = 1; i <= nclone; i++) {
        sum[i] += R[i];
        /* My (Susanta) adjustment to the problem of "sum" getting negative in some elements  */
        if (sum[i] <= 0) {
            if (i == 1) sum[i] = 1;
            else sum[i] = sum[i - 1];
        }
        //if(sum[i] <= 0) sum[i] = sum[i-1];
        ret -= log(sum[i]);
        //fprintf(myout,"R[i] %lf sum %lf ret %lf\n",R[i],sum[i],ret);
        //fprintf(stdout,"R[i] %lf sum %lf ret %lf clone number %d \n",R[i],sum[i],ret,i);
    }
    //fclose(myout);
    ret += const1;
    Free_Dec(sum, 1);

    return ret;

}

/*
 * Returns the gradient of the function
 */
void Gradient(double *gradient, double *y, double **aa, double *R,
        int nprobe, int nclone, int M, int *order) {

    double grad, den, num;
    int l, i, j;
    double *den_vec;

    den_vec = Dec_Vector(1, nclone);
    for (i = 1; i <= nclone; i++) {
        den_vec[i] = 0;
        for (j = 1; j <= nprobe + 1; j++)
            den_vec[i] -= (1 - aa[i][order[j]])*(1 - aa[i][order[j - 1]]) * Min(y[j], (double) M);
    }

    for (l = 1; l <= nprobe; l++) {
        grad = 0.0;
        for (i = 1; i <= nclone; i++) {
            den = den_vec[i];
            den += R[i];

            num = -Indicator(y[l], (double) M)*(aa[i][order[l]] - 1)*(aa[i][order[l - 1]] - 1) -
                    (aa[i][order[nprobe]] - 1) * Indicator(y[nprobe + 1], (double) M);
            /* A fix by Susanta: if den == 0.0 then fix den at den  = 0.001  */
            if (den == 0.0) den = 0.001;
            grad += num / den;
        }
        gradient[l] = -grad;
    }

    Free_Dec(den_vec, 1);



}

/*
 * Function to project the direction vector whenever it is necessary
 * in order to satisfy the constraints of the problem
 */
void Project(double *direction, double *y, int M, int nprobe) {
    int i, nzero, *flag, bool1;
    double last_old_direction;


    flag = Int_Vector(1, nprobe);
    for (i = 1; i <= nprobe; i++)
        flag[i] = 0;
    nzero = 0;

    do {
        bool1 = 1;
        last_old_direction = direction[nprobe + 1];
        for (i = 1; i <= nprobe; i++) {
            if ((y[i] < 0.0001) && (direction[i] < 0)) {
                bool1 = 0;
                direction[nprobe + 1] += direction[i];
                direction[i] = 0;
                flag[i] = 1;
                nzero++;
            } else if ((fabs(y[i]-(double) M) < 0.0001) && (direction[i] > 0)) {
                bool1 = 0;
                direction[nprobe + 1] += direction[i];
                direction[i] = 0;
                flag[i] = 1;
                nzero++;
            }

        }

        if (((y[nprobe + 1] < 0.0001) && (direction[nprobe + 1] < 0)) ||
                ((fabs(y[nprobe + 1]-(double) M) < 0.0001) && (direction[nprobe + 1] > 0))) {
            bool1 = 0;
            direction[nprobe + 1] /= (nprobe - nzero);
            for (i = 1; i <= nprobe; i++) {
                if (direction[i] != 0 && flag[i] != 1) direction[i] += direction[nprobe + 1];

            }
            direction[nprobe + 1] = 0;
        }
    } while (bool1 == 0);

    Free_Int(flag, 1);

}

/*
 * Function to reverse the current ordering for the probes
 */
void Reverse(int *piorder, int nprobe) {
    int n1, n2, temp, nn, j;

    do {
        // n1 = 1+(int)(nprobe*drand48()); USING int rand(void): range 0 to RAND_MAX
        //n1 = 1+(int)(nprobe*(double) (rand()/RAND_MAX));
        n1 = 1 + (int) (nprobe * ((double) rand() / RAND_MAX));
        //n2 = 1+(int)((nprobe-1)*(double) (rand()/RAND_MAX));
        n2 = 1 + (int) ((nprobe - 1)*(double) ((double) rand() / RAND_MAX));
        if (n2 >= n1) ++n2;
        else {
            temp = n1;
            n1 = n2;
            n2 = temp;
        }
        nn = n2 - n1 + 1;
    } while (nn == nprobe);
    nn = nn / 2;
    for (j = 1; j <= nn; j++) {
        temp = piorder[n1 + j - 1];
        piorder[n1 + j - 1] = piorder[n2 + 1 - j];
        piorder[n2 + 1 - j] = temp;
    }

}

void Print_Results(int *order, double *y, int nprobe, double func) {

    int i;


    fprintf(OUT, "\nThe final value of the objective function = %lf\n", func);
    fprintf(OUT, "The probe ordering is :\n\n");
    for (i = 1; i <= nprobe; i++)
        fprintf(OUT, "%d\n", order[i]);
    fprintf(OUT, "\n\nThe spacings are :\n\n");
    for (i = 1; i <= nprobe + 1; i++)
        fprintf(OUT, "%lf\n", y[i]);

    fclose(OUT);

}
